 // ==+==+==+==+==+==+==+==+==+==+==+==+==//
 // Modelo de Simulacion Sardina  (MOSAR) 
 // VERSION 1.0     2011-10-25
 //  LA Cubillos - EPOMAR
 // Universidad de Concepcion
 // Demandante: Subsecretaria de Pesca
 //======================================
 //Bitacora LOG
 //

DATA_SECTION
	!! cout<<"MOSAR detecta la plataforma "<<PLATFORM<< endl;

	init_adstring DataFile;
	init_adstring ControlFile;

	!! BaseFileName=stripExtension(ControlFile);
	!! cout<<BaseFileName<<endl;
	!! ReportFileName = BaseFileName + adstring(".rep");

	!! ad_comm::change_datafile_name(DataFile);

	int SimFlag;
	int rseed;
	LOC_CALCS
		SimFlag=0;
		rseed=999;
		int on,opt;
		//La linea siguiente evalua si la opcion de linea de comando "-SimFlag" existe
		//Si esta existe usa el numero aleatorio para la simulacion
		if((on=option_match(ad_comm::argc,ad_comm::argv,"-sim",opt))>-1)
		{
			SimFlag=1;
			rseed=atoi(ad_comm::argv[on+1]);
			cout<<"======================\n"<<endl;
			cout<<"      ** Simulacion Poblacional **"<<endl;
			cout<<"======================"<<endl;
		}
	END_CALCS

	//Lee los objetos del archivo de datos desde initprefix
	//Debug
	init_int debug
	//Time indicators
	init_int nmo  //number of months (time index)
	init_int nyr //number of years
	init_int styr //start year
	init_int endyr //end year
	!! cout<<"styr\t"<<styr<<endl;
	!! cout<<"endyr\t"<<endyr<<endl;
	!! cout<<"nmo\t"<<nmo<<endl;

	//Read time-series data from Tindex = 1 to nmo
	init_matrix s_TSdata(1,nmo,1,11);
	//indices and dependent
	vector Tindex(1,nmo);   //Indicador de tiempo secuencial
	vector Month(1,nmo);    //Meses
	vector Years(1,nmo);     //anhos
	//vector Effort(1,nmo);
	vector s_obs_catch(1,nmo);   //captura de sardina
	vector s_cpue(1,nmo);
	vector s_obs_Brec(1,nmo); //RECLAS acoustic biomass survey
	vector s_obs_Bpel(1,nmo); //PELACES acoustic biomass survey
	vector s_obs_SSB(1,nmo);  //DEPM spawning biomass
	vector s_aRLP(1,nmo);    //seasonal length-weight parameters
	vector s_bRLP(1,nmo);
	vector s_tempMat(1,nmo);  //seasonal maturity

        //Traspasa datos a vectores 
	LOC_CALCS
		//sardina
		Tindex = column(s_TSdata,1);  //Time index
		Years =  column(s_TSdata,2);   //Years
		Month = column(s_TSdata,3);   //Months
		//Effort = column(s_TSdata,8);
		s_obs_catch = column(s_TSdata,4);   //total monthly fishery catch biomass
		s_cpue = column(s_TSdata,5);    //Catch per unit effort
		s_obs_Brec=column(s_TSdata,6);  //Obs acoustic biomass Reclas survey
		s_obs_Bpel=column(s_TSdata,7); //Obs acoustic biomass Pelaces survey
		s_obs_SSB=column(s_TSdata,8);  //DEPM biomass
		s_aRLP=column(s_TSdata,10);     //intercept length-weight relationship
		s_bRLP=column(s_TSdata,11);     //slope LWR
		s_tempMat=column(s_TSdata,9);   //seasonal maturity
	END_CALCS

	//Length indicators
 	init_number s_stlen    //first length class
 	init_number s_steplen  //size step for length classes
 	init_number s_nlen     //number of length classes

	//Read fishery length composition
	init_int s_nobs_fishlen;
	init_ivector s_mo_fishlen(1,s_nobs_fishlen);
	init_matrix s_obs_p_fishlen(1,s_nobs_fishlen,1,s_nlen);
	//init_matrix s_obs_p_fishlen(1,nmo,1,s_nlen)

	!!if (debug==1){char mark; cout<<"s_catch= "<<s_obs_catch<<endl; cin>>mark; }  
	!!if (debug==1){char mark; cout<<"s_cpue= "<<s_cpue<<endl; cin>>mark; }  

	//Read Reclas length composition
	init_int s_nobs_reclaslen;
	init_ivector s_mo_reclaslen(1,s_nobs_reclaslen);
	init_matrix s_obs_p_reclaslen(1,s_nobs_reclaslen,1,s_nlen)

	//Read Pelaces length composition
	init_int s_nobs_pelacelen;
	init_ivector s_mo_pelacelen(1,s_nobs_pelacelen);
	init_matrix s_obs_p_pelacelen(1,s_nobs_pelacelen,1,s_nlen)

	vector s_len(1,s_nlen) //marca de clase de longitud

	//BIOLOGY
	init_int s_nages; //number of ages (yearly)
  	init_number s_M;        //Natural mortality rate annual
  	init_number s_Linf;
  	init_number s_Kappa;
  	init_number s_t0;

  	//init_number len1; //First length at age 1
  	init_number age1;				//recruitment age (3 mo for october, 4 mo for november...
  	init_vector s_CVageParms(1,2);	     //Intra-cohort parameters for CV of length at age
  	init_vector s_mat(1,2);	      		     //maturity parameters
  	//init_vector s_fec(1,2);     		     //batch fecundity parameters
	init_int nco; 				     	//number of cohorts per year during the time window
	!! cout<<"NCohorts\t"<<nco<<endl;
	init_ivector mes_rec(1,nco);		//initial month of recruitment cohorts
	init_ivector mes_surv(1,nco);		//last month of surviving cohorts
	init_ivector minageco(1,nco);		//minima edad por cohortes

	!!if (debug==1){char mark; cout<<"minageco= "<<minageco<<endl; cin>>mark; }  

	//some additional
	matrix ages(1,12,1,s_nages);
	number s_len1;
	int j;
	int i;
	int k;
	int l;
	number pi;
	int itmp;
	int surv; //index of survival month of cohort
	int rec;  //index of initial month of cohorts
	imatrix mes(1,12,1,nyr)
	matrix s_ageco(1,nmo,1,nco);  //ages per cohort
	matrix s_lenco(1,nmo,1,nco);   //length at age per cohort
	matrix s_weico(1,nmo,1,nco);	    //weight at age per cohort
	matrix s_matco(1,nmo,1,nco);  //maturity at age per cohort
	matrix ones(1,nmo,1,nco);  //units
	vector tspawn(1,nyr);      //month peak spawning
	vector treclas(1,nyr);	   //month reclas survey
	vector tpelace(1,nyr);    //month pelaces
        vector tfebre(1,nyr);
	vector tmarzo(1,nyr); //mes en marzo
	vector s_w_age(1,s_nages)
	vector s_mat_age(1,s_nages)

        number s_rho;
	number s_rho_m
	matrix s_sd_age(1,12,1,s_nages)
	matrix s_cv_age(1,12,1,s_nages)
	number diff
	3darray s_Agesize(1,12,1,s_nages,1,s_nlen)
	matrix s_len_age(1,12,1,s_nages)
	number s_Mm; //monthly natural mortality
	number s_nagesm

	LOC_CALCS
		s_Mm=s_M/12;    //Fixed nat mort per month
		s_nagesm=s_nages*12;  //ages months
		//Compute mean length for length-classes
		s_len(1) = s_stlen;
		for(j=1;j<s_nlen;j++)
		{
			s_len(j+1)=s_len(j)+s_steplen;
		}

		//time index for peak spawning biomass (=mes 11=agosto because 1=october)
		//to be used in stock-recruitment
		tspawn(1)=11;	//August
		treclas(1)=4;  //January
		tpelace(1)=8; //May
		tfebre(1)=5;  //febrero
		tmarzo(1)=6; //marzo
		for(i=1;i<nyr;i++)
		{
			tspawn(i+1)=tspawn(i)+12;
			treclas(i+1)=treclas(i)+12;
			tpelace(i+1)=tpelace(i)+12;
			tfebre(i+1)=tfebre(i)+12;
			tmarzo(i+1)=tmarzo(i)+12;
		}
	END_CALCS
	!!if (debug==1){char mark; cout<<"slen"<<s_len<<endl; cin>>mark; }  

	!!if (debug==1){char mark; cout<<"tspawn"<<tspawn<<endl; cin>>mark; }  


	LOC_CALCS
		//age, length-at-age, w-at-age, maturity at age
		for(i=1; i<=nco;i++)
		{
			surv=mes_surv(i);
			rec=mes_rec(i);
			itmp=minageco(i);
			for(j=rec;j<=surv;j++)
			{
				s_ageco(j,i) = itmp;	
				s_lenco(j,i) = s_Linf*(1-exp(-1*s_Kappa*(s_ageco(j,i)-s_t0)/12));
				s_weico(j,i) = s_aRLP(j)*pow(s_lenco(j,i),s_bRLP(j));
				s_matco(j,i) = 1./(1+exp(-(s_mat(1)+s_mat(2)*s_lenco(j,i))));
				ones(j,i)=1.;
				itmp=itmp+1;
			}
		}
	END_CALCS


	LOC_CALCS
		//AGE LENGTH STUFF
		s_len1=s_Linf*(1-exp(-1*s_Kappa*(age1-s_t0)));
		s_len_age(1,1)=s_len1;
		ages(1,1)=age1;
		s_rho=mfexp(-1.*s_Kappa);
		//ages
		//First month
		for(j=1;j<s_nages;j++)
			{
					s_len_age(1,j+1)=s_Linf*(1.-s_rho)+s_rho*s_len_age(1,j);
					ages(1,j+1)=ages(1,j)+1;
			}
		//Now for the rest 11 months
		s_rho_m=mfexp(-1.*s_Kappa/12);
		for(j=1;j<=s_nages;j++)
			{
			for(i=1;i<12;i++)
				 {
					s_len_age(i+1,j)=s_Linf*(1.-s_rho_m)+s_rho_m*s_len_age(i,j);
					ages(i+1,j)=ages(i,j)+0.0833;
				}
			}
			//CV of length at age
			for(i=1;i<=12;i++)
				{
				for(j=1;j<=s_nages;j++)
					{
						s_cv_age(i,j)=(s_CVageParms(1)+s_CVageParms(2)*ages(i,j))/100;
						if(s_cv_age(i,j) < 0){s_cv_age(i,j)=0.002324; }
						//s_cv_age(i,j)=0.05;
						s_sd_age(i,j)=s_cv_age(i,j)*s_len_age(i,j);
					}
				}

			//calcula avg-wt-age and maturity for august peak spw s-r stuff age-structured
			for (j=1;j<=s_nages;j++)
			{
				s_w_age(j) = s_aRLP(11)*pow(s_len_age(11,j),s_bRLP(11));
				s_mat_age(j)=s_tempMat(11)/(1+exp(-(s_mat(1)+s_mat(2)*s_len_age(11,j))));
       			}

			//AGE-LENGTH KEY
			     pi = 3.141592653590;
			     for(k=1;k<=12;k++)
				{
				for(i=1;i<=s_nages;i++)
						{
						for(j=1;j<=s_nlen;j++)
							{
								diff = s_len(j)-s_len_age(k,i);
							     s_Agesize(k,i,j)=(1/sqrt(2*pi*(square(s_sd_age(k,i)))))*mfexp(-0.5*square(diff)/(square(s_sd_age(k,i))));
							}
							s_Agesize(k,i)=s_Agesize(k,i)/sum(s_Agesize(k,i));
						}
				}

				for(i=1;i<=nyr;i++)
				{
					for(j=1;j<=12;j++)
					{
						mes(j,i)=(i-1)*12+j;
					}
				}

					 // cout << mes << endl;
	END_CALCS

	!!if (debug==1){char mark; cout<<"watage= "<<s_w_age<<endl; cin>>mark; }  
	!!if (debug==1){char mark; cout<<"matage= "<<s_mat_age<<endl; cin>>mark; }  

	//End of data file
	init_int eof;	
	LOC_CALCS
	  cout<<"eof = "<<eof<<endl;
	  if(eof==999){
		cout<<"\n ** FIN DE LA SECCION DATA ** \n"<<endl;
	  }else{
		cout<<"\n ¡¡¡¡ ERROR AL LEER LOS DATOS!!!!! \n"<<endl; exit(1);
	  }
	END_CALCS

	number fmsy; //Fishing mortality at Fmsy
	number msy; //Maximum sustainable yield
	number bmsy; //Spawning biomass at MSY
	vector size_tau(1,3); //MLE of variance for size composition
	matrix simcatsize(1,nmo,1,s_nlen); //captura por talla para la simulacion

	//============================================
	// Lee los parametros de control de estimacion
	//============================================
	//Leading parameters
	//theta[1]		log_ro, or log_msy     unexploited recruitment
	//theta[2]		steepness(h), or log_fmsy
	//theta[3]		log_m  : mortality rate
	//theta[4]		log_avgrec : average recruitment
	//theta[5]		prho  : proportion of observation variance recruitment
	//theta[6]		varphi :precision

	!! ad_comm::change_datafile_name(ControlFile);
	init_int npar;
	init_matrix theta_control(1,npar,1,7)
	
	!!if (debug==1){char mark; cout<<"Control= "<<theta_control<<endl; cin>>mark; }  

	vector theta_ival(1,npar);
	vector theta_lb(1,npar);
	vector theta_ub(1,npar);
	ivector theta_phz(1,npar);
	ivector theta_prior(1,npar);

	number th_ival;
	number th_lb;
	number th_ub;
	number th_phz;
	number th_prior;

	LOCAL_CALCS
		theta_ival = column(theta_control,1);  //ival 
		theta_lb = column(theta_control,2);
		theta_ub = column(theta_control,3);
		theta_phz = ivector(column(theta_control,4)); //phase
		theta_prior = ivector(column(theta_control,5)); //prior

	END_CALCS

	init_int ngear; //numero de arte con length frequency
	//=========================
	//       Opciones de selectividad
	//=========================
	// type 1 = logistic (2 parms)
	// type 2 = logistica fija y predeterminada cambiando la phase (sel_phz) a negativa
	init_ivector isel_type(1,ngear);   //tipo de selectividad para la pesqueria, reclas y pelaces
	init_vector ahat(1,ngear);           //age-at-50% vulnerability logistic
	init_vector ghat(1,ngear);	       //std at 50%
	init_ivector sel_phz(1,ngear);	//Phase for estimating selectivity parameters.
	!! cout<<isel_type<<endl;
        ivector isel_npar(1,ngear);

	// fixed logistic (no parameters estimated)
	// ensure sel_phz is set to negative value.
	LOCAL_CALCS
		isel_npar.initialize();
		for(i=1;i<=ngear;i++)
		{
			switch(isel_type(i))
			{
				case 1:
				     //logistic
				     isel_npar(i)=2;
				     break;
				case 2:
				     //logistic fixed
				     isel_npar(i)=2;
				     if(sel_phz(i)>0) sel_phz(i) = -1;
				     break;
				case 3:
				     //just another
				default:break;
			}

		}
		cout<<"Number of estimated selectivity parameters\n"<<isel_npar<<endl;
	END_CALCS

	//Controls for prior on catchability q (including cpue).
	init_int nits;       //numero de surveys (inits = 4, 1: cpue; 2: reclas; 3: pelaces; 4: mph)
	init_ivector q_prior(1,nits);  //tipo de prior usado para q(0: uniform; 1: normal;..)
	init_vector q_mu(1,nits);   //valor de q assumido escala log
	init_vector q_sd(1,nits);   // std para q
	//!! cout<<"q Prior\n"<<q_mu<<endl<<q_sd<<endl;

	init_vector cv_surveys(1,nits); //CVs surveys
	init_vector nsizemult(1,ngear) //effective size for multinomial
	number CV_cpue;
	number CV_reclas;
	number CV_pelaces;
	number CV_mph;
	number n_fish;
	number n_reclas;
	number n_pelaces
	LOCAL_CALCS
		CV_cpue = cv_surveys(1);
		CV_reclas = cv_surveys(2);
		CV_pelaces = cv_surveys(3);
		CV_mph = cv_surveys(4);
		n_fish = nsizemult(1);
		n_reclas = nsizemult(2);
		n_pelaces = nsizemult(3);
	END_CALCS

	//Miscellaneous controls
	// 1 -> verbose
	// 2 -> recruitment model (1=beverton-holt, 2=rickers)
	// 3 -> autocorrelation in recruitment
	// 4 -> std in catch first phase
	// 5 -> std in catch in last phase
	// 6 -> mean fishing mortality rate to regularize the solution
	// 7 -> standard deviation of mean F penalty in first phases
	// 8 -> standard deviation of mean F penalty in last phase.

	
	init_vector cntrl(1,8);
	int verbose;

	init_int eofc;
	LOC_CALCS
		verbose = cntrl(1);
		cout<<"cntrl\n"<<cntrl<<endl;
		cout<<"eofc\t"<<eofc<<endl;
		if(eofc==999){
			  cout<<"\n ** FIN CONTROL ** \n"<<endl;
	  	}else{
			  cout<<"\n ¡¡¡¡ ERROR CONTROL!!!!! \n"<<endl; exit(1);
	 	}
	END_CALCS

	int nf;
	
	ivector ilvec(1,6);
	!!ilvec=ngear;			//number of gear for selectivity ngear
	!!ilvec(2)=nits;		//number of surveys antes nit
	!!ilvec(3)= ngear;         //na_gears number of size-comps
	!!ilvec(4)=1;

PARAMETER_SECTION

	init_bounded_number_vector theta(1,npar,theta_lb,theta_ub,theta_phz);
	!! for(int i=1;i<=npar;i++) theta(i)=theta_ival(i);




	!!if (debug==1){char mark; cout<<"theta= "<<theta<<endl; cin>>mark; }  

	//catchability
	init_number_vector log_q(1,nits,3)
	//!!for(i=1;i<=nits;i++) log_q(i) = q_mu(i);

	!!if (debug==1){char mark; cout<<"log_q= "<<log_q<<endl; cin>>mark; }

	//selectivity
	init_bounded_number_vector sel_par(1,ngear,-5,5,sel_phz); 
	init_bounded_number_vector sel_ds(1,ngear,-5,5,sel_phz); 

	LOCAL_CALCS
		for(k=1;k<=ngear;k++)
		{
			if(isel_type(k)==1 || isel_type(k)==2)
			{
				sel_par(k) = log(ahat(k));
				sel_ds(k) = log(ghat(k));
			}
		}
	END_CALCS
	!!if (debug==1){char mark; cout<<"sel_par= "<<sel_par<<endl; cin>>mark; }
	!!if (debug==1){char mark; cout<<"sel_ds= "<<sel_ds<<endl; cin>>mark; }

	//Fishing mortality rate parameters
	init_bounded_vector log_ft_pars(1,nmo,-30.,3.0,1);

	LOC_CALCS
		if(!SimFlag) log_ft_pars = log(0.045);
	END_CALCS
	//Annual recruitment deviations
	init_bounded_vector log_rec_devs(1,nco,-15.,15.,2);

	objective_function_value f;

	!!if (debug==1){char mark; cout<<"Obj fun= "<<f<<endl; cin>>mark; }

	//DEPENDENCIAS
	number ro;		//unfished recruits
	number bo;		//unfished spawning stock biomass
	number kapsr;
	number m;		//initial natural mortality rate
	number log_avgrec;	//log of average recruitment.
	number prho;		//proportion of the observation error
	number varphi		//total precision in the CPUE & Rec anomalies.
	number so;
	number beta;
	vector ssb(1,nyr);		//spawning stock biomass yearly
	vector rt(2,nyr); 		//predicted  recruits from S-R curve yearly
	vector delta(2,nyr);	//residuals for annual stock recruitment yearly
	vector log_rt(1,nco); //log recruitment

	vector q(1,nits);			//survey catchability coefficients
	matrix sel_fish(1,nmo,1,nco);      //selectivity fishery
	matrix sel_reclas(1,nmo,1,nco);       //selectivity survey reclas
	matrix sel_pelaces(1,nmo,1,nco);   //selectivity survey pelaces
	vector avg_log_sel(1,3); //conditional penalty for objective function

       //mortality rate
	vector ft(1,nmo);
	matrix F(1,nmo,1,nco);	  //Fishing mortality per cohorts
	matrix M(1,nmo,1,nco);
	matrix Z(1,nmo,1,nco);
	matrix S(1,nmo,1,nco);

	//population
	matrix N(1,nmo,1,nco);     //Number per cohorts
	matrix sb(1,nmo,1,nco);  //spawning biomass
	matrix biom(1,nmo,1,nco);
        matrix npop(1,nmo,1,s_nages);
	matrix npopsize(1,nmo,1,s_nlen);
        vector ntot(1,nmo); //total population
	vector btot(1,nmo);
        vector sbiom(1,nmo);
	vector vulnbiom(1,nmo); //vulnerable biomass

	//ingresar predic	
	matrix s_catco(1,nmo,1,nco); //predicted catch-per-cohort
	matrix s_pred_p_co(1,nmo,1,nco); //predicted catch prop per cohort
	vector s_pred_catch(1,nmo); //predicted monthly catch
	matrix s_catage(1,nmo,1,s_nages);
	matrix s_pred_p_age(1,nmo,1,s_nages);
	matrix s_catlen(1,nmo,1,s_nlen);
	matrix s_pred_p_fishlen(1,nmo,1,s_nlen);
	vector s_pred_cpue(1,nmo)
	vector s_pred_Brec(1,nmo)
	vector s_pred_Bpel(1,nmo)
	vector s_pred_SSB(1,nmo)
	matrix s_pred_p_reclasage(1,nmo,1,s_nages)
	matrix s_pred_p_pelaceage(1,nmo,1,s_nages)
	matrix s_pred_p_reclaslen(1,nmo,1,s_nlen)
	matrix s_pred_p_pelacelen(1,nmo,1,s_nlen)

	vector offset(1,ngear);
	matrix nlvec(1,6,1,ilvec);	//matrix for negative loglikelihoods

	!!if (debug==1){char mark; cout<<"nlvec= "<<nlvec<<endl; cin>>mark; }


PRELIMINARY_CALCS_SECTION

  //Calcula offset para el multinomial
  for (int i=1; i<= s_nobs_fishlen; i++)
  {
	s_obs_p_fishlen(i)=s_obs_p_fishlen(i)/sum(s_obs_p_fishlen(i));
	for (int j=1;j <= s_nlen; j++)
	{
		if (s_obs_p_fishlen(i,j)>0.0)
		{
			offset(1)-=n_fish*s_obs_p_fishlen(i,j)*log(s_obs_p_fishlen(i,j));
       		}
	}
   }

  for (int i=1; i<= s_nobs_reclaslen; i++)
  {
	s_obs_p_reclaslen(i)=s_obs_p_reclaslen(i)/sum(s_obs_p_reclaslen(i));
	for (int j=1;j <= s_nlen; j++)
	{
		if (s_obs_p_reclaslen(i,j)>0.0)
		{
			offset(2)-=n_reclas*s_obs_p_reclaslen(i,j)*log(s_obs_p_reclaslen(i,j));
       		}
	}
   }

  for (int i=1; i<= s_nobs_pelacelen; i++)
  {
	s_obs_p_pelacelen(i)=s_obs_p_pelacelen(i)/sum(s_obs_p_pelacelen(i));
	for (int j=1;j <= s_nlen; j++)
	{
		if (s_obs_p_pelacelen(i,j)>0.0)
		{
			offset(3)-=n_reclas*s_obs_p_pelacelen(i,j)*log(s_obs_p_pelacelen(i,j));
       		}
	}
   }

	//Corre el modelo con parametros conocidos para simular datos
	nf=0;
	if(SimFlag)
	 {
		initParameters();
		calcSelectivities();
		calcTotalMortality();
		//simulation_model(rseed);
	  }

	//if (debug==1){char mark; cout<<"DFLpelace= "<<s_obs_p_pelacelen<<endl; cin>>mark; }


RUNTIME_SECTION
    maximum_function_evaluations 100,200,500,25000,25000
    convergence_criteria 0.01,0.01,1.e-5,1.e-5

PROCEDURE_SECTION
	
	initParameters();
	calcSelectivities();
	calcTotalMortality();
	calcNumbersAtAge();
        calcFisheryObservations();		
	calc_survey_observations();
	
	calc_stock_recruitment();
	
	calc_objective_function();

        //exit(1);

	//sd_depletion=sbt(nyr)/bo;
	
	//if(mceval_phase()) mcmc_output();


FUNCTION initParameters
  {
	/*
	Esta funcion permite extraer los parametros para ser usados.
	
	Esta rutina debe correr antes del modelo de simulacion.
	*/
	
	ro = mfexp(theta(1));
	dvariable h = theta(2);
	switch(int(cntrl(2)))
	{
		case 1:
			//Beverton-Holt model
			kapsr = (4.*h/(1.-h));
			break;
		case 2:
			//Ricker model
			kapsr = pow((5.*h),1.25);
		break;
	}
		
	//TODO Alternative parameterization using MSY and FMSY as leading parameters
	//TODO autocorrelation in recruitment
	
	m = mfexp(theta(3));
	log_avgrec = theta(4);
	
	prho=theta(5);
	varphi=theta(6);
	if(verbose)cout<<"**** Ok after initParameters ****"<<endl;
  }

FUNCTION calcSelectivities
  {
	/*
		This function loops over each ngear and calculates the corresponding
		selectivity for that gear type. It first uses a switch statement 
		to calculate selectivity curves based on isel_type where:
		1) logistic selectivity with 2 parameters
		2) Fixed logistic
	*/
	int k,i,j;
	dvariable p1,p2;
	sel_fish.initialize();
	sel_reclas.initialize();
	sel_pelaces.initialize();
	
	for(k=1;k<=ngear;k++) //antes ngear
	{
		switch(isel_type(k))
		{
			case 1:
				// logistic selectivity for case 1 or 2
				p1 = mfexp(sel_par(k));
				p2 = mfexp(sel_ds(k));
				if(k==1){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_fish(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}

				if(k==2){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_reclas(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}
				if(k==3){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_pelaces(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}

				break;
			
			case 2:
				p1 = mfexp(sel_par(k));
				p2 = mfexp(sel_ds(k));
				if(k==1){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_fish(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}

				if(k==2){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_reclas(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}
				if(k==3){
					for(i=1;i<=nco;i++)
					{
						rec=mes_rec(i);
						surv=mes_surv(i);
						for(j=rec;j<=surv;j++)
						{
							sel_pelaces(j,i)=1./(1.+mfexp((p1-s_ageco(j,i))/p2));
						}
					}
				}

				break;
				
			case 3:
				// age-specific selectivity coefficients
				break;
				
			//default:
				//break;
				
		}  // switch
	
	}

	if(verbose)cout<<"**** Ok after calcSelectivities ****"<<endl;	
  }

FUNCTION calcTotalMortality
  {
	/*
	This routine calculates fishing mortality, total mortality
	and annaul survival rates (exp(-Z)) for each cohort in each
	month.
	
	There is a complication in that if there is a survey gear
	then the fishing mortality rate is set to an extremely small
	value: exp(-70.)~3.975e-31, which is effectively 0.
	
	The above issue is no longer an issue b/c now estimating Fs
	for years where there is non-zero catch.
	*/

	int i,j;
	F.initialize();
	ft.initialize();
	
	//Fishing mortality
	for(j=1;j<=nmo;j++)
 	{
		ft(j)=mfexp(log_ft_pars(j));
		for (int i=1;i<=nco;i++)
		{
			F(j,i)=sel_fish(j,i)*ft(j);
       		}
	}
		
	//Natural mortality
	//M(1,nmo,1,nco);
	M = ones*m;
	M = M/12;
	Z=M+F;
	S=mfexp(-Z);


	if(verbose) cout<<"**** OK after calcTotalMortality ****"<<endl;
	
  }

FUNCTION calcNumbersAtAge
  {
     int i,j;
     N.initialize();

     //Estimate recruitment
     log_rt = log_avgrec+log_rec_devs;
     for(i=1;i<=nco;i++)
	{
		rec=mes_rec(i);
		N(rec,i)=mfexp(log_rt(i));
	}
	//Now the entire population matrix
     for(i=1;i<=nco;i++)
	{
		rec=mes_rec(i);
		surv=mes_surv(i);
		for(j=rec;j<surv;j++)
		{
			N(j+1,i)=N(j,i)*S(j,i);
		}
	}

     for(i=1;i<=nco;i++)
	{
		for(j=1;j<=nmo;j++)
		{
			biom(j,i)=N(j,i)*s_weico(j,i);   //biomass per cohort
			sb(j,i)=biom(j,i)*s_matco(j,i); //mature biomass per cohort
		}
	}
     for(j=1;j<=nmo;j++)
	{
		ntot(j)=0;
		btot(j)=0;
		sbiom(j)=0;
		vulnbiom(j)=0;	
		for(i=1;i<=nco;i++)
		{
			ntot(j)+=N(j,i);
			btot(j)+=biom(j,i);
			sbiom(j)+=sb(j,i);
			vulnbiom(j)+=sel_fish(j,i)*biom(j,i)*mfexp(-0.5*Z(j,i));
		}
		sbiom(j)=sbiom(j)*s_tempMat(j); //mature biomass corrected by seasonal maturity
	}

        //if (debug==1){char mark; cout<<"Sbiom= "<<sbiom<<endl; cin>>mark; }

       if(verbose)cout<<"**** Ok after calcNumbersAtAge ****"<<endl;
  }

FUNCTION calcFisheryObservations
  {
	int j,i,itmp;
	for(i=1;i<=nco;i++)
	{
		rec=mes_rec(i);
		surv=mes_surv(i);
		for(j=rec;j<=surv;j++)
		{
			s_catco(j,i)=F(j,i)*N(j,i)*(1.-S(j,i))/Z(j,i); //catch per cohort
		}
	}
	
	for(i=1;i<=nmo;i++)
	{
		s_pred_p_co(i)=s_catco(i)/sum(s_catco(i));
		s_pred_catch(i)=0;
		itmp=0;
		for(j=1;j<=nco;j++)
		{
			s_pred_catch(i)+=s_catco(i,j)*s_weico(i,j);		
			//ordena catco en catage
			if(ones(i,j)>0)
			{
				itmp=itmp+1;
				s_catage(i,itmp)=s_catco(i,j);
				s_pred_p_age(i,itmp)=s_pred_p_co(i,j); //
				npop(i,itmp)=N(i,j);
			}
		}
	}


	 //Predicted catch and population at length
	 for (i=1;i<=nyr;i++)
	 {
		for(j=1;j<=12;j++)
		{
			itmp=mes(j,i);
			s_catlen(itmp)=s_catage(itmp)*s_Agesize(j);
			npopsize(itmp)=npop(itmp)*s_Agesize(j);
      		}
	}
	for (i=1; i<= nmo; i++)
	{
		s_pred_p_fishlen(i)=s_catlen(i)/(sum(s_catlen(i))+1e-10);
	}

        //if (debug==1){char mark; cout<<"Catage= "<<s_catage<<endl; cin>>mark; }

        if(verbose)cout<<"**** Ok after calcFisheryObservations ****"<<endl;
  }

FUNCTION calc_survey_observations
  {
	int i,j,k,ii;
	//qs: 1= cpue; 2: reclas; 3: pelaces; 4: mph
	q=mfexp(log_q);
	//pred_cpue
	s_pred_cpue=q(1)*vulnbiom;

	//Reclas survey and pelaces
	for(i=1;i<=nmo;i++)
	{
		s_pred_Brec(i)=0;
		s_pred_Bpel(i)=0;
		for(j=1;j<=nco;j++)
		{
			s_pred_Brec(i)+=q(2)*sel_reclas(i,j)*s_weico(i,j)*N(i,j);
			s_pred_Bpel(i)+=q(3)*sel_pelaces(i,j)*s_weico(i,j)*N(i,j);
		}
    	}
	 //DEPM ssb
	 s_pred_SSB=q(4)*sbiom;

	 //prepare to compute survey at length
	for(i=1;i<=nmo;i++)
	{
		itmp=0;
		for(j=1;j<=nco;j++)
		{
			if(ones(i,j)>0)
			{
				itmp=itmp+1;
				s_pred_p_reclasage(i,itmp)=sel_reclas(i,j)*N(i,j);
				s_pred_p_pelaceage(i,itmp)=sel_pelaces(i,j)*N(i,j); //
			}
		}
				s_pred_p_reclasage(i)= s_pred_p_reclasage(i)/(sum(s_pred_p_reclasage(i))+1e-10);
				s_pred_p_pelaceage(i)= s_pred_p_pelaceage(i)/(sum(s_pred_p_pelaceage(i))+1e-10);
	}

	 //length comp surveys 2 and 3
	 for (i=1;i<=nyr;i++)
	 {
		for(j=1;j<=12;j++)
		{
				itmp=mes(j,i);
				s_pred_p_reclaslen(itmp)=s_pred_p_reclasage(itmp)*s_Agesize(j);
				s_pred_p_pelacelen(itmp)= s_pred_p_pelaceage(itmp)*s_Agesize(j);
		}
	}

        //if (debug==1){char mark; cout<<"qs= "<<q<<endl; cin>>mark; }

	if(verbose)cout<<"**** Ok after calc_survey_observations ****"<<endl;

  }

FUNCTION calc_stock_recruitment
  {
	/*
	The following code is used to derive unfished
	spawning stock biomass bo and the stock-
	recruitment parameters for the:
	Beverton-Holt Model
		-Rt=k*Ro*St/(Bo+(k-1)*St)*exp(delta-0.5*tau*tau)
	Ricker Model
		-Rt=so*St*exp(-beta*St)*exp(delta-0.5*tau*tau)
		
	
	CHANGED Need to adjust spawning biomass to post fishery numbers.
	CHANGED Need to adjust spawners per recruit (phib) to average fecundity.
	*/ 
	int i,itmp;
	dvariable tau = (1.-prho)/varphi;
	dvar_vector tmp_rt(2,nyr);
	dvar_vector lx(1,s_nages); lx=1;
        for(i=1+1;i<=s_nages;i++) lx(i)=lx(i-1)*exp(-m);

        dvariable phib;
	for(i=1;i<=s_nages;i++)
        {
          phib+= lx(i)*s_w_age(i)*s_mat_age(i);	//
	}
        // if (debug==1){char mark; cout<<"phib= "<<phib<<endl; cin>>mark; }

	dvariable so = kapsr/phib;		//max recruits per spawner
	dvariable beta;
	bo = ro*phib;  					//unfished spawning biomass
         //if (debug==1){char mark; cout<<"bo= "<<bo<<endl; cin>>mark; }
	
	//sbt=rowsum(elem_prod(N,fec));			//SM Dec 6, 2010
	//CHANGED adjusted spawning biomass downward by ctrl(12)
	//AQUI HAY QUE CAMBIAR Y TRAER LA BIOM DESOV POR AÑO desde sbiom y usando peak spaw
	//OJO REVISAR tspawn y que hay delay en un year
	for(i=1;i<=nyr;i++)
	{
		itmp = tspawn(i);
		ssb(i)=sbiom(itmp);
	}
         if (debug==1){char mark; cout<<"ssb= "<<ssb<<endl; cin>>mark; }

	
	dvar_vector tmp_st(1,nyr-1);
	for(i=1;i<nyr;i++)
	{
		tmp_st(i) = ssb(i);
	}


	switch(int(cntrl(2)))
	{
		case 1:
			//Beverton-Holt model
			beta = (kapsr-1.)/bo;
			for(i=1;i<nyr;i++)
			{
				//tmp_rt=elem_div(so*tmp_st,1.+beta*tmp_st);
				tmp_rt(i+1)=so*tmp_st(i)/(1.+beta*tmp_st(i));
			}
			break;
		case 2:
			//Ricker model
			//dvariable so=kapsr/phib;
			beta=log(kapsr)/bo;
			for(i=1;i<nyr;i++)
			{
				//tmp_rt=elem_prod(so*tmp_st,exp(-beta*tmp_st));
				tmp_rt(i+1)=so*tmp_st(i)*mfexp(-beta*tmp_st(i));
			}

		break;
	}
         if (debug==1){char mark; cout<<"temprt= "<<tmp_rt<<endl; cin>>mark; }

	//residuals in stock-recruitment curve
        dvar_vector tmrt(2,nyr);
        for(i=2;i<=nyr;i++)
	{
          tmrt(i)=log_rt(i+s_nages-1);
  	}

	rt=exp(tmrt); //trans(N)(1)(syr+1,nyr);
	delta = log(rt)-log(tmp_rt)+0.5*tau*tau; //residuals from S-R

       if (debug==1){char mark; cout<<"logrt= "<<log_rt<<endl; cin>>mark; }

        if (debug==1){char mark; cout<<"rt= "<<rt<<endl; cin>>mark; }

       if (debug==1){char mark; cout<<"delta= "<<delta<<endl; cin>>mark; }
	
	if(verbose)cout<<"**** Ok after calc_stock_recruitment ****"<<endl;
	
  }

FUNCTION calc_objective_function
  {
	/*q_prior is an ivector with current options of 0 & 1.
	0 is a uniform density (ignored) and 1 is a normal
	prior density applied to log(q).*/

	/*
	There are several components to the objective function
	Likelihoods:
		-1) likelihood of the catch data
		-2) likelihood of the survey abundance index
		-3) likelihood of the survey size comps
		-4) likelihood of the fishery size comps
		-5) likelihood for stock-recruitment relationship
		-6) likelihood for fishery selectivities
	*/
        int i,j;
        int ii;
        nlvec.initialize();
        //1) likelihood of the catch data
	double sig_c =cntrl(4);
	if(last_phase()) sig_c=cntrl(5);
	if(active(log_ft_pars))
	for(i=1;i<=nmo;i++)
	{
		if(s_obs_catch(i)!=0)
		nlvec(1,1)+=dnorm(log(s_pred_catch(i)),log(s_obs_catch(i)),sig_c);
	}

       if (debug==1){char mark; cout<<"nlvec catch= "<<nlvec(1,1)<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"sigC= "<<s_obs_catch<<endl; cin>>mark; }

        //2) likelihood of the survey data
        //cpue
        for(i=1;i<=nmo;i++)
	{
		if(s_cpue(i)!=0)
		nlvec(2,1)+=dnorm(log(s_pred_cpue(i)),log(s_cpue(i)),CV_cpue);
	}
       if (debug==1){char mark; cout<<"nlvec cpue= "<<nlvec(2,1)<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"obs cpue= "<<s_cpue<<endl; cin>>mark; }

       //reclas
        for(i=1;i<=s_nobs_reclaslen;i++)
	{
		ii=s_mo_reclaslen(i);
		nlvec(2,2)+=dnorm(log(s_pred_Brec(ii)),log(s_obs_Brec(ii)),CV_reclas);
	}
       if (debug==1){char mark; cout<<"nlvec reclas= "<<nlvec(2,2)<<endl; cin>>mark; }

       //pelaces
        for(i=1;i<=s_nobs_pelacelen;i++)
	{
		ii=s_mo_pelacelen(i);
		nlvec(2,3)+=dnorm(log(s_pred_Bpel(ii)),log(s_obs_Bpel(ii)),CV_pelaces);
	}
       if (debug==1){char mark; cout<<"nlvec pelaces= "<<nlvec(2,3)<<endl; cin>>mark; }

       for(i=1;i<=nmo;i++)
	{
		if(s_obs_SSB(i)!=0)
		nlvec(2,4)+=dnorm(log(s_pred_SSB(i)),log(s_obs_SSB(i)),CV_mph);
	}
       if (debug==1){char mark; cout<<"nlvec mph= "<<nlvec(2,4)<<endl; cin>>mark; }

        //3) and 4) likelihood of the fishery and survey size comp data

       //if (active(rec_dev))
        //{
              //rec_like=lambda(1)*norm2(rec_dev);
              nlvec(3,1)=0.;
              //f+=rec_like;
              for (i=1; i<=s_nobs_fishlen; i++)
              {
                       ii=s_mo_fishlen(i);
                       for (j=1; j<=s_nlen; j++)
                       {
                         nlvec(3,1)-=n_fish*s_obs_p_fishlen(i,j)*log(s_pred_p_fishlen(ii,j)+1.e-13);
                       }
              }
              nlvec(3,1)-=offset(1);

              nlvec(3,2)=0.;
              for (i=1; i<=s_nobs_reclaslen; i++)
              {
                       ii=s_mo_reclaslen(i);
                       for (j=1; j<=s_nlen; j++)
                       {
                         nlvec(3,2)-=n_reclas*s_obs_p_reclaslen(i,j)*log(s_pred_p_reclaslen(ii,j)+1.e-13);
                       }
              }
              nlvec(3,2)-=offset(2);

              nlvec(3,3)=0.;
              for (i=1; i<=s_nobs_pelacelen; i++)
              {
                       ii=s_mo_pelacelen(i);
                       for (j=1; j<=s_nlen; j++)
                       {
                         nlvec(3,3)-=n_pelaces*s_obs_p_pelacelen(i,j)*log(s_pred_p_pelacelen(ii,j)+1.e-13);
                       }
              }
              nlvec(3,3)-=offset(3);


        //}
       if (debug==1){char mark; cout<<"nlvec fishsize= "<<nlvec(3,1)<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"nlvec reclassize= "<<nlvec(3,2)<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"nlvec pelacesize= "<<nlvec(3,3)<<endl; cin>>mark; }

	//5) likelihood for stock-recruitment relationship
	dvariable tau = (1.-prho)/varphi;
	if(active(theta(1)))
		nlvec(4,1)=dnorm(delta,tau);

       if (debug==1){char mark; cout<<"prho= "<<prho<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"tau= "<<tau<<endl; cin>>mark; }
       if (debug==1){char mark; cout<<"nlvec= "<<nlvec(4,1)<<endl; cin>>mark; }

	/*
	PRIORS for estimated model parameters from the control file.
	*/
	dvar_vector priors(1,npar);
	dvariable ptmp; priors.initialize();
	for(i=1;i<=npar;i++)
	{
		if(active(theta(i)))
		{
			switch(theta_prior(i))
			{
			case 1:		//normal
				ptmp=dnorm(theta(i),theta_control(i,6),theta_control(i,7));
				break;
				
			case 2:		//lognormal CHANGED RF found an error in dlnorm prior. rev 116
				ptmp=dlnorm(theta(i),theta_control(i,6),theta_control(i,7));
				break;
				
			case 3:		//beta distribution (0-1 scale)
				double lb,ub;
				lb=theta_lb(i);
				ub=theta_ub(i);
				ptmp=dbeta((theta(i)-lb)/(ub-lb),theta_control(i,6),theta_control(i,7));
				break;
				
			case 4:		//gamma distribution
				ptmp=dgamma(theta(i),theta_control(i,6),theta_control(i,7));
				break;
				
			default:	//uniform density
				ptmp=1./(theta_control(i,3)-theta_control(i,2));
				break;
			}
			priors(i)=ptmp;	
		}
	}




	//Priors for suvey q based on control file.
	dvar_vector qvec(1,nits);
	qvec.initialize();
	for(i=1;i<=nits;i++)
	{
		if(q_prior(i)==2) 
		{
			qvec(i)=dnorm(log(q(i)),q_mu(i),q_sd(i));
		}
	}

       if (debug==1){char mark; cout<<"qvec= "<<qvec<<endl; cin>>mark; }

	/*
	The following are penalties that are invoked in early
	phases to ensure reasonable solutions are obtained,
	and to reduce the sensitivity of initial starting
	conditions.  Penalties include:
		-1) keep average fishing mortality rate near 
			0.2 and in the last phase relax this constraint.
		-3) normal prior for log rec devs with std=50.
	*/
	
	dvar_vector pvec(1,2);
	pvec.initialize();
	
	//Penalties to regularize the solution for fishing mortality rates
	dvariable log_fbar = mean(log_ft_pars);
	if(last_phase())
	{
		pvec(1) = dnorm(log_fbar,log(cntrl(6)),cntrl(8));
		//Penalty for log_rec_devs (large variance here)
		pvec(2) = dnorm(log_rec_devs,2.0);
	}
	else
	{
		pvec(1) = dnorm(log_fbar,log(cntrl(6)),cntrl(7));
		//Penalty for log_rec_devs (CV ~ 0.0707) in early phases
		pvec(2)=100.*norm2(log_rec_devs);
	}


       if (debug==1){char mark; cout<<"pvec= "<<pvec<<endl; cin>>mark; }




      f=sum(nlvec)+sum(pvec)+sum(qvec); //+sum(priors);
      nf++;
  //end of likelihoods
    }




TOP_OF_MAIN_SECTION
	time(&start);
	arrmblsize = 50000000;
	gradient_structure::set_GRADSTACK_BUFFER_SIZE(1.e7);
	gradient_structure::set_CMPDIF_BUFFER_SIZE(1.e7);
	gradient_structure::set_MAX_NVAR_OFFSET(5000);
	gradient_structure::set_NUM_DEPENDENT_VARIABLES(5000);



GLOBALS_SECTION
	#undef REPORT
	#define REPORT(object) report << #object "\n" << object << endl;
    
	#undef COUT
	#define COUT(object) cout << #object "\n" << object <<endl;

	#if defined(WIN32) && !defined(__linux__)
		const char* PLATFORM = "Windows";
	#else
		const char* PLATFORM = "Linux";
	#endif

	#include <admodel.h>
	#include <time.h>
	#include <string.h>
	#include <stats.cxx>
	time_t start,finish;
	long hour,minute,second;
	double elapsed_time;

	 
	adstring BaseFileName;
	adstring ReportFileName;
	
	adstring stripExtension(adstring fileName)
	{
		/*
		This function strips the file extension
		from the fileName argument and returns
		the file name without the extension.
		*/
		const int length = fileName.size();
		for (int i=length; i>=0; --i)
		{
			if (fileName(i)=='.')
			{
				return fileName(1,i-1);
			}
		}
		return fileName;
	}

